###Step 1 Building New Dataset from Raw Data
nn <- 1 #number of file
Data <- paste("Old_",nn,".csv",sep='')
Data <- read.csv(Data)
head(Data)
Data1 <- subset(Data, subset = Filter == 1)
name <- Data$Name[1]

##Building Data and find n
dat <- Data1
dat[,6] <- dat[,6]+1
Data2 <- rbind(Data1,dat)
head(Data2)
tail(Data2)

#Using naturalSpline to fit the curve
#Library
library(splines2)

Phase <- Data2$Phase
Light <- Data2$Light

########3
n <- 2
repeat{
	sp_phase <- naturalSpline(Phase,intercept = TRUE, df = n)
	model <- lm(Light~ sp_phase)
	SumModel <- summary(model)
	Rsquare <- SumModel$r.squared
	if(Rsquare > 0.9999){
		break
	}
	n <- n+1
}
n

sp_phase <- naturalSpline(Phase,intercept = TRUE, df = n)
model <- lm(Light~ sp_phase)
SumModel <- summary(model)
SumModel

plot(Phase,Light, pch = 16, cex=0.5, ylim = c(0,1.5))
lines(Phase,model$fitted.values, type = 'l')

Coefficients <- as.numeric(model$coefficients)
	#sp_phase
	#FittedValues <- model$fitted.values
	#head(sp_phase)
	#head(FittedValues)

#Make new data from the curve
X <- seq(0,2,0.01)
Ycap <- predict(sp_phase, X)
m <- length(X)

NewPD <- c()
for(j in 1:m){
	Y <- as.numeric(Ycap[j,])
	Sum <- Coefficients[1]
	for(i in 2:(n-1)){
		a <- Coefficients[i]*Y[i-1]
		Sum <- Sum + a
	}
	NewPD <- c(NewPD,Sum)
}

plot(Phase,Light, pch = 16, cex=0.5)
lines(Phase,model$fitted.values, type = 'l')
lines(X,NewPD, col = 'red')
title(name)

Name <- rep(name,m)
NewData <- cbind(Name,X,NewPD)
namefile <- paste("New_",nn,".csv", sep="")
write.table(NewData,namefile, sep=",", col.names = NA)

nameplot <- paste(name)
savePlot(filename = nameplot, type = "bmp")

### Step 2 Finging W1, W2, D1 and D2 from New Dataset
Data <- paste("New_",nn,".csv",sep='')
Data <- read.csv(Data)
head(Data)

X <- Data$X
Y <- Data$NewPD
name <- Data[1,2]

#Finding Slope (DY)
n <- length(X)
id <- 1:n
DY <- NA
for(i in 2:n){
	s <- (Y[i]-Y[i-1])/(X[i]-X[i-1])
	DY <- c(DY,s)
}

#Finding difference of Slope :Diff.DY
Diff.DY <- c(NA,NA)
for(i in 3:n){
	d <- abs(DY[i]-DY[i-1])
	Diff.DY <- c(Diff.DY,d)
}

##Finding difference of Diff.DY :Diff2.DY
Diff2.DY <- c(NA,NA,NA)
for(i in 4:n){
	d <- abs(Diff.DY[i]-Diff.DY[i-1])
	Diff2.DY <- c(Diff2.DY,d)
}

NewData <- data.frame(X,Y,DY,Diff.DY,Diff2.DY,id)

namefile <- paste(name,"_new2.csv", sep="")
write.table(NewData,namefile, sep=",", col.names = NA)

###########################################################

############idxDY by DY(Slope +/-) and idxDY by DDY(+/-):DYC ## slope change point
idxDY <- c()
XDY <- c()
start <- sign(DY[2])
for(i in 3:n){
	if(sign(DY[i]) != start){
		idxDY <- c(idxDY,i)
		XDY <- c(XDY,X[i])
		start <- sign(DY[i])
	}
}

DataDYC <- NewData[idxDY,]

## SelectPoints
SelectPoints <- data.frame()
###########################################################
###########Find inflaction Points : IFP
P1 <- NewData[0:25,]
P2 <- NewData[26:50,]
P3 <- NewData[51:75,]
P4 <- NewData[76:100,]
P5 <- NewData[101:125,]
P6 <- NewData[126:150,]
P7 <- NewData[151:175,]
P8 <- NewData[176:201,]

IFP <- data.frame()

P1 <- P1[order(P1$DY),1:6]
n <- length(P1[,1])
IFP <- rbind(IFP,P1[n-1,])
P2 <- P2[order(P2$DY),1:6]
n <- length(P1[,2])
IFP <- rbind(IFP,P2[1,])

P3 <- P3[order(P3$DY),1:6]
n <- length(P3[,1])
IFP <- rbind(IFP,P3[n,])
P4 <- P4[order(P4$DY),1:6]
n <- length(P4[,2])
IFP <- rbind(IFP,P4[1,])

P5 <- P5[order(P5$DY),1:6]
n <- length(P5[,1])
IFP <- rbind(IFP,P5[n,])
P6 <- P6[order(P6$DY),1:6]
n <- length(P6[,2])
IFP <- rbind(IFP,P6[1,])

P7 <- P7[order(P7$DY),1:6]
n <- length(P7[,1])
IFP <- rbind(IFP,P7[n,])
P8 <- P8[order(P8$DY),1:6]
n <- length(P8[,2])
IFP <- rbind(IFP,P8[1,])

IFP
idx.IFP <- IFP$id
m <- length(IFP[,1])
NamePoint <- rep('IFP',m)
SelectPoints <- data.frame(NamePoint,IFP)

#############
idx.IFP

k <- length(idxDY)

Int1 <- idx.IFP[1]:idx.IFP[2]
Int2 <- idx.IFP[3]:idx.IFP[4]
Int3 <- idx.IFP[5]:idx.IFP[6]
Int4 <- idx.IFP[7]:idx.IFP[8]


### Vertex Points
DataVTP <- data.frame()
###1
m <- length(Int1)
VT1 <- c()
for(i in 1:m){
	for(j in 1:k){
		if(Int1[i] == idxDY[j]){
			VT1 <- c(VT1,Int1[i])
		}
	}
}
Data.p1 <- NewData[VT1,]
Data.p1 <- Data.p1[order(Data.p1$Y),1:6]
m <- length(Data.p1[,1])
DataVTP <- rbind(DataVTP,Data.p1[m,])

###2
m <- length(Int2)
VT2 <- c()
for(i in 1:m){
	for(j in 1:k){
		if(Int2[i] == idxDY[j]){
			VT2 <- c(VT2,Int2[i])
		}
	}
}
Data.p2 <- NewData[VT2,]
Data.p2 <- Data.p2[order(Data.p2$Y),1:6]
m <- length(Data.p2[,1])
DataVTP <- rbind(DataVTP,Data.p2[m,])

###3
m <- length(Int3)
VT3 <- c()
for(i in 1:m){
	for(j in 1:k){
		if(Int3[i] == idxDY[j]){
			VT3 <- c(VT3,Int3[i])
		}
	}
}
Data.p3 <- NewData[VT3,]
Data.p3 <- Data.p3[order(Data.p3$Y),1:6]
m <- length(Data.p3[,1])
DataVTP <- rbind(DataVTP,Data.p3[m,])

###4
m <- length(Int4)
VT4 <- c()
for(i in 1:m){
	for(j in 1:k){
		if(Int4[i] == idxDY[j]){
			VT4 <- c(VT4,Int4[i])
		}
	}
}
Data.p4 <- NewData[VT4,]
Data.p4 <- Data.p4[order(Data.p4$Y),1:6]
m <- length(Data.p4[,1])
DataVTP <- rbind(DataVTP,Data.p4[m,])

m <- length(DataVTP[,1])
NamePoint <- rep('VTP',m)
DataVTP <- data.frame(NamePoint,DataVTP)
SelectPoints <- rbind(SelectPoints,DataVTP)



#################The most still point :MSP
idx.IFP
idx.VTP <- DataVTP$id
idx.VTP

NP1 <- (idx.IFP[1]+1):(idx.VTP[1]-1)
NP2 <- (idx.VTP[1]+1):(idx.IFP[2]-1)

NP3 <- (idx.IFP[3]+1):(idx.VTP[2]-1)
NP4 <- (idx.VTP[2]+1):(idx.IFP[4]-1)

NP5 <- (idx.IFP[5]+1):(idx.VTP[3]-1)
NP6 <- (idx.VTP[3]+1):(idx.IFP[6]-1)

NP7 <- (idx.IFP[7]+1):(idx.VTP[4]-1)
NP8 <- (idx.VTP[4]+1):(idx.IFP[8]-1)

MSP <- data.frame()
##NP1
DataNP1 <- NewData[NP1,]
DataNP1 <- DataNP1[order(DataNP1$Diff2.DY),1:6]
MSP <- rbind(MSP,DataNP1[1,])

##NP2
DataNP2 <- NewData[NP2,]
DataNP2 <- DataNP2[order(DataNP2$Diff2.DY),1:6]
MSP <- rbind(MSP,DataNP2[1,])

##NP3
DataNP3 <- NewData[NP3,]
DataNP3 <- DataNP3[order(DataNP3$Diff2.DY),1:6]
MSP <- rbind(MSP,DataNP3[1,])

##NP4
DataNP4 <- NewData[NP4,]
DataNP4 <- DataNP4[order(DataNP4$Diff2.DY),1:6]
MSP <- rbind(MSP,DataNP4[1,])

##NP5
DataNP5 <- NewData[NP5,]
DataNP5 <- DataNP5[order(DataNP5$Diff2.DY),1:6]
MSP <- rbind(MSP,DataNP5[1,])

##NP6
DataNP6 <- NewData[NP6,]
DataNP6 <- DataNP6[order(DataNP6$Diff2.DY),1:6]
MSP <- rbind(MSP,DataNP6[1,])

##NP7
DataNP7 <- NewData[NP7,]
DataNP7 <- DataNP7[order(DataNP7$Diff2.DY),1:6]
MSP <- rbind(MSP,DataNP7[1,])

##NP8
DataNP8 <- NewData[NP8,]
DataNP8 <- DataNP8[order(DataNP8$Diff2.DY),1:6]
MSP <- rbind(MSP,DataNP8[1,])

m <- length(MSP[,1])
NamePoint <- rep('MSP',m)
MSP <- data.frame(NamePoint,MSP)
SelectPoints <- rbind(SelectPoints,MSP)


#############################################################
###Find Mean Cutoff Points
Mean <- mean(c(MSP$Y,DataVTP$Y))

MCOP <- c()

#Interval 1 VTP1 & IFP2
D1 <- NewData[SelectPoints$id[9]:SelectPoints$id[2],]
DiffFromMean <- abs(D1$Y-Mean)
D1 <- data.frame(D1,DiffFromMean)
D1 <- D1[order(D1$DiffFromMean),1:6]
MCOP <- rbind(MCOP,D1[1,])

#Interval 2 IFP3 & VTP2
D2 <- NewData[SelectPoints$id[3]:SelectPoints$id[10],]
DiffFromMean <- abs(D2$Y-Mean)
D2 <- data.frame(D2,DiffFromMean)
D2 <- D2[order(D2$DiffFromMean),1:6]
MCOP <- rbind(MCOP,D2[1,])

#Interval 3 VTP2 & IFP4
D3 <- NewData[SelectPoints$id[10]:SelectPoints$id[4],]
DiffFromMean <- abs(D3$Y-Mean)
D3 <- data.frame(D3,DiffFromMean)
D3 <- D3[order(D3$DiffFromMean),1:6]
MCOP <- rbind(MCOP,D3[1,])

#Interval 4 IFP5 & VTP3
D4 <- NewData[SelectPoints$id[5]:SelectPoints$id[11],]
DiffFromMean <- abs(D4$Y-Mean)
D4 <- data.frame(D4,DiffFromMean)
D4 <- D4[order(D4$DiffFromMean),1:6]
MCOP <- rbind(MCOP,D4[1,])

m <- length(MCOP[,1])
NamePoint <- rep('MCOP',m)
MCOP <- data.frame(NamePoint,MCOP)
SelectPoints <- rbind(SelectPoints,MCOP)

############### Find Lowest Points: LWP
LWP <- c()
#Interval_1 IFP2 & IFP3
D1 <- NewData[SelectPoints$id[2]:SelectPoints$id[3],]
D1 <- D1[order(D1$Y),1:6]
LWP <- rbind(LWP,D1[1,])

#Interval_2 IFP4 & IFP5
D2 <- NewData[SelectPoints$id[4]:SelectPoints$id[5],]
D2 <- D2[order(D2$Y),1:6]
LWP <- rbind(LWP,D2[1,])

m <- length(LWP[,1])
NamePoint <- rep('LWP',m)
LWP <- data.frame(NamePoint,LWP)
SelectPoints <- rbind(SelectPoints,LWP)

###############
m <- length(SelectPoints[,1])
Name <- rep(name,m)
SelectPoints <- cbind(Name,SelectPoints)

namefile <- paste(nn,"_SelectPoints.csv", sep="")
write.table(SelectPoints ,namefile, sep=",", col.names = NA)

#############################################################
##Final Plot
plot(X,Y,pch = 16,cex = 0.3, ylim = c(0,1.5), main = name)
lines(X,Y, col = 'red')
abline(h = Mean, col = 'green')

abline(v= IFP$X)
abline(v = MSP$X, lty = 2, col = 'brown')
abline(v= DataVTP$X, lty=3, col= 'dark green')
points(MSP$X,MSP$Y, pch=16,cex=0.7)
points(DataVTP$X,DataVTP$Y, pch=16,cex=0.7)

points(X,Diff.DY,pch = 16,cex = 0.3, col = 'blue')
points(X,Diff2.DY,pch = 16,cex = 0.3, col = 'brown')

nameplot <- paste(name)
savePlot(filename = nameplot, type = "bmp")

########################################
#####Collect data
 Name <- c()
 D1.ver <- c()
 D2.ver <- c()
 W1.ver <- c()
 W2.ver <- c()
###########################################
for(nn in 1:50){

Data <- paste(nn,"_SelectPoints.csv",sep = '')
Data <- read.csv(Data)
Data
Data1 <- Data[9:20,]
Mean <- mean(Data1$Y)

d1 <- Mean-Data[25,5]
d2 <- Mean-Data[26,5]
w1 <- Data[22,4]-Data[21,4]
w2 <- Data[24,4]-Data[23,4]
if(d1 > d2){
	D1 <- d1
	D2 <- d2
	W1 <- w1
	W2 <- w2
}else{
	D1 <- d2
	D2 <- d1
	W1 <- w2
	W2 <- w1
}
name <- Data[1,2]
###########################################

Name <- c(Name,name)
D1.ver <- c(D1.ver,D1)
D2.ver <- c(D2.ver,D2)
W1.ver <- c(W1.ver,W1)
W2.ver <- c(W2.ver,W2)

}
############################################
final <- data.frame(Name,D1.ver,D2.ver,W1.ver,W2.ver)
final
namefile <- paste("Final_3.csv", sep="")
write.table(final,namefile, sep=",", col.names = NA)

###Step 3 Model A
for(k in 1:5){

#filter
f <- k
name <- paste("Final_Data F",f,".csv",sep = "")
data <- read.csv(name)
ratio.W <- data$W1/data$W2
ratio.D <- data$D1/data$D2
ratio.T <- data$T1/data$T2
ratio.Om <- data$Om1/data$Om2
Data <- data.frame(data[,1:6],ratio.W,ratio.D,data[7:22],ratio.T,ratio.Om)
head(Data)

### Multiple Regression

Coef <- rbind()
Res <- c()
SIG <- c()

ResponVar <- c()
Rsquare <- c()
F.Stat <- c()
Pvalue <- c()
Significant <- c()

D1 <- Data$D1
D2 <- Data$D2
W1 <- Data$W1
W2 <- Data$W2

for(i in 9:26){
	responVar <- names(Data)[i]
	model <- lm(Data[,i] ~ D1+D2+W1+W2)
	sm <- summary(model)

	res <- c()
	Sig <- c()
	for(j in 1:5){
		res <- c(res,responVar)
		if(sm$coefficients[, "Pr(>|t|)"][j] < 0.05){
			Sig <- c(Sig,"sig")
		}else{
			Sig <- c(Sig,"not sig")
		}
	}
	
	#18*5
	Coef <- rbind(Coef,sm$coefficients)
	Res <- c(Res,res)
	SIG <- c(SIG,Sig)

	fstat <- sm$fstatistic
	pval_F <- pf(fstat[1], fstat[2], fstat[3], lower.tail = FALSE)

	#18	
	ResponVar <- c(ResponVar,responVar)
	Rsquare <- c(Rsquare,sm$r.squared)
	F.Stat <- c(F.Stat,sm$fstatistic[1])
	Pvalue <- c(Pvalue,pval_F)

	if(pval_F < 0.05){
		Significant <- c(Significant,"sig")
	}else{
		Significant <- c(Significant,"not sig")
	}
}

OP1 <- data.frame(Res,Coef,SIG)
OP1
OP2 <- data.frame(ResponVar,Rsquare,F.Stat,Pvalue,Significant )
OP2

n1 <- paste("Output1 A_F",f,".csv",sep='')
n2 <- paste("Output2 A_F",f,".csv",sep='')
write.csv(OP1,n1)
write.csv(OP2,n2)

}

###Step 4 Model B
for(k in 1:5){

#filter
f <- k
name <- paste("Final_Data F",f,".csv",sep = "")
data <- read.csv(name)
ratio.W <- data$W1/data$W2
ratio.D <- data$D1/data$D2
ratio.T <- data$T1/data$T2
ratio.Om <- data$Om1/data$Om2
Data <- data.frame(data[,1:6],ratio.W,ratio.D,data[7:22],ratio.T,ratio.Om)
head(Data)

### Multiple Regression

Coef <- rbind()
Res <- c()
SIG <- c()

ResponVar <- c()
Rsquare <- c()
F.Stat <- c()
Pvalue <- c()
Significant <- c()

D <- ratio.D
W <- ratio.W

for(i in 9:26){
	responVar <- names(Data)[i]
	model <- lm(Data[,i] ~ D + W)
	sm <- summary(model)

	res <- c()
	Sig <- c()
	for(j in 1:3){
		res <- c(res,responVar)
		if(sm$coefficients[, "Pr(>|t|)"][j] < 0.05){
			Sig <- c(Sig,"sig")
		}else{
			Sig <- c(Sig,"not sig")
		}
	}
	
	#18*5
	Coef <- rbind(Coef,sm$coefficients)
	Res <- c(Res,res)
	SIG <- c(SIG,Sig)

	fstat <- sm$fstatistic
	pval_F <- pf(fstat[1], fstat[2], fstat[3], lower.tail = FALSE)

	#18	
	ResponVar <- c(ResponVar,responVar)
	Rsquare <- c(Rsquare,sm$r.squared)
	F.Stat <- c(F.Stat,sm$fstatistic[1])
	Pvalue <- c(Pvalue,pval_F)

	if(pval_F < 0.05){
		Significant <- c(Significant,"sig")
	}else{
		Significant <- c(Significant,"not sig")
	}
}

OP1 <- data.frame(Res,Coef,SIG)
OP1
OP2 <- data.frame(ResponVar,Rsquare,F.Stat,Pvalue,Significant )
OP2

n1 <- paste("Output1 B_F",f,".csv",sep='')
n2 <- paste("Output2 B_F",f,".csv",sep='')
write.csv(OP1,n1)
write.csv(OP2,n2)

}
